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Abstract 


Component adaptive grid (CAG) methods for solving hyperbolic partial differen- 
tial equations (PDE’s) are discussed in this paper. Applying recent stability results 
for a class of numerical methods on uniform grids, the convergence of these methods 
for linear problems on component adaptive grids is established here. Furthermore, 
the computational error can be estimated on CAGs using the stability results. Using 
these estimates, the error can be controlled on CAG’s. Thus, the solution can be 
computed efficiently on CAG’s within a given error tolerance. Computational results 
for time dependent linear problems in one and two space dimensions are presented. 
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1 Introduction 


Component adaptive grid methods for solving hyperbolic PDE’s were introduced in 
the early 1980’s. An overview of the method is given in Section 2. More details can be 
found in Berger and Oliger [2], and Berger [1]. However, the grid structure used in this 
paper is different from the one Berger used. As discussed in Section 2, stair step grids 
like those of Chesshire and Henshaw [4] are used in our CAG methods here, instead 
of rotated rectangular grids. One major component of the adaptive strategy is to 
estimate the local truncation error at each grid point, then refine where the estimated 
errors are larger than a given tolerance 8. The smaller 8 is, the smaller the final error 
e is expected to be in some weighted L 2 norm. However, no quantitative relationship 
between these two kinds of errors had been established. Recently, new stability results 
have been developed by Pelle Olsson [6] which allow us to establish such a relationship 
for large classes of problems and methods. The results can be applied to various classes 
of problems, e.g., those of hyperbolic, parabolic and hyperbolic-parabolic type, using 
a large class of numerical methods on uniform grids. As we will see in Section 3, the 
structures of component adaptive grids allow us to define the solution on piecewise 
uniform grids. So the stability theories can be applied on CAG’s. Convergence for 
linear problems using these methods on CAG’s is proved in Section 3. Also the 
tolerance 8 on local truncation error is estimated in term of the tolerance e on the 
final error. Furthermore, the results in Section 3 will also help us estimate the final 
error using simple quadrature, and serve us as guidelines on developing strategies 
for CAG methods, since we have a very good understanding of the sources and the 
magnitudes of various computational errors. Finally, some computional results for 
time dependent problems in one and two space dimensions are given in Section 4. 


2 An Overview of Component Adaptive Grids 

We first introduce some notation for our discussion. Suppose the problem we wish to 
solve is written as 


U t = 

Lu + f 

on ft x [0, T] 

(1) 

u(0) = 

u 0 

on ft 

(2) 

B u = 

b 

on dQ, x [0, T] 

(3) 


where ft C R d is a bounded domain in physical space, L is a spatial partial differetial 
operator on ft and tiGif 1 . We assume this to be a well-posed initial-boundary value 
problem which is defined in Section 3. Let ft^, dCth and [0, T]k be the discretizations of 
ft, dft and [0, T], respectively. In Section 3, these discretizations are defined precisely 
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for our component adaptive grids. For the time being, we can consider them as 
general grids. 

Let u/, be a grid function defined on f l h x [0, T]*. We will discuss the use of 
finite difference methods on these grids. Without loss of generality, and avoiding 
complicated notation, we write our methods in explicit one-step form as 


v h (t - 1 - k) 

= L h v h {t) + k f h (t) 

on n h x [o,r] fc 

(4) 

MO) 

— U'Oh, 

on f lh. 

(5) 

Bh Vh(i ) 

= b h (t) 

on dClh. x [0, T]k 

(6) 


where we use subscripts to denote projections of functions onto the appropriate grids 
and discretizations of operators on these grids. If Uh is the projection of the exact 
solution of the above system onto f lh, then 

u h (t + k) = L h u h {t) + kf h (t) + kr h on f)/, x [0, T) k (7) 

where r* is the local truncation error. This notation will also be used on the piecewise 
uniform grids which we will discuss next. 


2.1 Composite Grids 

In real applications, the physical domains often have complicated geometries. In order 
to use finite difference schemes on these domains, we decompose the physical domain 
and transform the parts into computational domains. However this topic is not the 
focus of this paper. Here only a brief introduction is given to make our presentation 
self contained. Details can be found in Chesshire and Henshaw [4], Venkata, Oliger 
and Ferziger [8], and Venkata [9]. 

We begin by forming a base composite grid 


G , o = UG°,j (8) 

3 

which will be characterized by a discretization parameter h 0 . 

This is well illustrated in Figure 1 where Go consists of the component grids Go.i, 
G 0 , 2 , and Go, 3. Go, 2 is a stair step grid with grid lines parallel to the coordinate axes. 
Such grids are called regular grids. 


Definition 1: A regular grid is a connected stair step grid of uniformly spaced 

points in each coordinate direction, and its grid lines are parallel to the coordinate 
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axes in either physical or computational space. 

No coordinate transformation is needed to solve the equation on regular grids in 
physical space. The curvilinear grids Go,i and Go ,3 are defined by specifying their 
boundaries and cuts. Regular girds in computational space are then mapped onto 
these grids in physical space using coordinate tranformations. To reduce clutter in 
Figure 1 , grids Go, 1 and Go , 3 are shown only in computational space. The component 
grids are chosen to obtain a sufficiently accuate representation of dQ, by di lh- ho is 
an estimate of the step size required to obtain a sufficiently accuate approximation of 
the solution over at least some specified fraction of the domain. The difficult problem 
here is to generate grids on the boundaries. The Bezier family of curves and surfaces 
are used to generate boundary grids in 2-D and 3-D, respectively (see [8] and [9]). 



Figure 1: Adaptive Composite Grid Structure 


2.2 Component Adaptive Grids 

The component grids mentioned in Section 2.1 are specified to describe the domain 
and its boundaries. Next, we will discuss the use of adaptive grids which are created 
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and destroyed during the course of the computation in order to maintain sufficient 
accuracy throughout the entire domain. 

During this process, we will create L— 1 additional refinement levels of composite 
grids G\ on top of the base grid Go- So the whole grid system can be written as 
follows: 

tj G ‘ (9) 

/=0 

where at each level /, we have 

G, = \JG U . (10) 

These will have spatial discretization parameters hi = ho/m l , where m might range 
from 2 to 10 depending on the problems being solved. According to Berger and Oliger 
[2], m = 4 is a reasonable choice for many hyperbolic problems. We usually let L be 
2, 3 or 4. As dictated by the nature of the problem and numerical algorithm, we main- 
tain an appropriate relationship between the spatial and the temporal discretization 
parameters, hi and ki. In particular, for problems which are essentially hyperbolic in 
character, we usually use the same mesh ratio on all levels, i.e., 

A = ki/hi = constant. (11) 

Another very important feature for our adaptive grids is that the grids are level 

nested, i.e., the region Gi is fully contained within the region G/_i. See the shaded 
refined grids G\,i, G 1,2, G\ ,3 and G\ ,4 in Figure 1. Since these grids have the same mesh 
ratio, they have the space-time structure illustrated in Figure 2. We take several time 
steps on the finer grid for each time step on the coarser grid. 



x 

Figure 2: Space-time grid structure 


We generate the refinement grids in response to computed estimates of local trun- 
cation error. If the solution is sufficiently smooth, a variant of Richardson extrapola- 
tion, call step- doubling, can be used to estimate the local truncation error at each 
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point of the grid. For simplicity, we consider the one-step explicit difference scheme 
described in equations (4) to (7). We define the operator Qh as follows: 


QhUh{x,t) = L h u h {x,t) + kf h {t). (12) 

If the solution is smooth enough, then the local truncation error can be written as 

kT K (x,t) = u h (x,t + k) - QhUh(x,t) 

= k(k qi a(x, t) + h q2 b(x, t)) + kO(k qi+1 + h q2+1 ) 

= Jfcr + kO(k qi+1 + h q2+1 ) (13) 

where the leading term is denoted by kr. If u is smooth enough and we take two time 
steps with the operator Qh, to leading order the error is 2 kr, 

Uh(x,t + 2k) — Q\uh{x,t) = 2 kr + kO(k qi+1 + h q2+1 ). (14) 

Let Q 2 h be the same difference operator as Qh but based on mesh widths of 2 h and 
2k. Also, assume the order of accuracy in time and space are equal, q 1 = q 2 = q. 
Then 


Uh(x,t + 2k) - Q 2h Uh{x,t) = 2k((2k) q a(x,t) + (2h) q b(x,t)) + 0(h q+2 ) 

= 2 q+1 kr + 0(h q+2 ). (15) 


Subtract equation (14) from equation (15) and use equation (13) to obtain 


n(x,t) 


Qlu h (x,t)-Q 7h u h {x,t) +1 

k( 2**' - 2 ) +0(/ “ > 


(16) 


The restriction that the accuracy in time and space is the same is not a severe one. 
For details see Berger [1]. For nonsmooth solutions we no longer have an accurate 
error estimate. However, the Richardson estimates still provide a good criterion for 
refinement since the estimates will be large near a singularity. 


2.3 Adaptive Grid Generation and Integration 

With the basic background mentioned above, we will now explain the adaptive grid 
generation and integation processes. It is important to organize the data structure 
for the adaptive composite grid in terms of connected components, i.e., connected 
sets of component grids at each level. See Figure 3 for the tree representing the grid 
structure in Figure 1. Details on implementation can be found in [1] and [2]. 
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Figure 3: Tree structure of connected components 


Evaluation of h 0 . Having begun with an estimated initial value of /io, we perform 
a trial integration on this grid and estimate the error. If more than a given 
fraction, say 1/2, of the grid fails to meet the error criterion and needs to be 
refined, we refine the whole grid with ho reduced by a factor of m, repeating this 
process if necessary. (This overall refinement is also done during the course of 
the computation if necessary.) The fraction 1/2 is used because of the overhead 
associated with grid refinement. It has been ascertained experimentally that 
about 3/4 of a domain can be refined (with the associated overhead) at the 
same computational cost as needed to compute on the entire refined domain 
without the adaptive method’s overhead. Once h 0 is established, the other hi 
are defined in terms of it. 

Time integration. The solution is advanced in time as follows. The basic operation 
is to solve on a connected component. A time step of k 0 is first taken on Go, 
and then each grid G/, l = 1, 2, . . . , L — 1 is in turn advanced by &/, connected 
component by connected component. Then Gi-\ is brought up to t + &l -2 and 
so on, until all of the grids are brought up to t + k 0 . Whenever a refinement is 
brought up to the time level of the next coarser grid, the point values on the 
coarser grid under the refinement are replaced by the solution on the refinement. 
Values of the solution on the interior edges of connected components at each 
level are obtained by interpolating in time the values already obtained on the 
next coarser grid. It should become clear at the end of Section 3 that this 
interpolation needs to have a certain order of accuracy to maintain the optimal 
rate of convergence. 

Grid modification. Local truncation errors are estimated by the process of step* 
doubling at fixed numbers of time steps on each grid level, usually every 4 to 
8 steps depending on the size of buffer zone mentioned below. Points at which 
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the error criterion is near violation are flagged, clusters of the flagged points are 
formed, and refinements overlaying the clusters are constructed. Once a new 
component grid is defined, the solution at previously refined points is retained, 
and the solution at all other points is obtained by interpolation accurate to a 
certain order from the parent grid. (At the initial time, in contrast, the solution 
at each refinement level is obtained from the initial data ti 0 .) The refined grids 
are constructed with buffer zones around the flagged points, sufficiently wide 
so that a phenomenon requiring refinement cannot move out of the refined area 
before the next regridding. This is necessary to maintain accuracy. The grids 
are moved by constructing a new grid and deleting the old one, see Figure 2. 
Each refinement level is in turn refined as necessary in the same manner, until 
the finest level has been reconstructed. 


3 Stability, Convergence and Error Estimation 

3.1 Stability 

As mentioned in the beginning of section 2, the original and discretized problems we 
deal with are formulated in equations (1) - (3) and equations (4) - (6), respectively. 
Once again, we recall that subscripts h and k are used to represent the discretized 
domains, operators and variables. Our goal is to use CAG’s to compute the solution 
such that the error is bounded by a given tolerance e, i.e., 

IMr)lk < ^ (u) 

where || • ||n h is some discrete norm and eh{t) is the discrete error function, both to 
be defined later. 

Before discussing the discrete case, we want to assume that the original initial- 
boundary problem is well-posed. The well-poseness of a large class of problems is 
discussed in Kreiss and Lorenz [5]. For hyperbolic and parabolic problems, we use 
the following definition. 


Definition 2: The initial-boundary value problem defined by equations (1) - (3) is 

well-posed if there exists an unique solution which satisfies the following estimate 

< Ke"‘(IM-)IIS + ll/(-, -JHLio., ] + ll i ’(-.')lllnx|o,ij) (18) 

where K and a are constants, t 6 [0, T] and the norms are defined in the usual way: 

IK-iOlln = / W(x,t)\ 2 dx (19) 

Jn 
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(20) 

( 21 ) 


ll&(v)llinx[o,t] = / I \b(x,r)\ 2 ds dr 

0 Oil 

ll/(') *)llnx[0,t] = J o f \f(x,T)\* dx dr 
and | • | is the standard L 2 norm in /2 n , i.e., \u\ 2 = u*u. 

In order to achieve our goal, some stability results, which are similar to the well- 
poseness of the original problem, are needed for the finite difference schemes. In his 
report, Pelle Olsson [6] used energy methods and summation by parts (a discrete 
version of integration by parts) to establish stability estimates on uniform grids for 
the semi-discretized case, i.e., when only the spatial domain is discretized. The spa- 
tial discretizations used are centered differences with various orders of accuracy. To 
overcome the difficulty of general boundary conditions, some projection matrices are 
used so that the solution can be projected into the solution space where the bound- 
ary conditions are satisfied. These results hold for a large class of linear problems, 
e.g., hyperbolic, parabolic and some mixed types. Although only 2-D problems were 
considered in his report, the results can be easily generalized to higher dimensions. 
We assume the following stability estimates holds for our finite difference methods on 
the uniform mesh. This is consistent with the results of Olsson [6]. 


Definition 3: A finite difference method for the initial-boundary value problem 

(equations (4) to (6)) on a uniform grid is stable if 

ll u fc(*)lln h < A'e at (||uo h ||n h + ||/fc||n hX [o,i] t + ll^llan h x[o,t] fc ) ( 22 ) 

where 1< and a are independant of k and h, t € [0, T] and the discrete norms are 
weighted L 2 norms in general. Here we can assume they have the following form: 

IM0II&, = E ( 23 > 

IWIkxM. = E kh "-' (24) 

IIAIIn.xM, = E (25) 

xjefWielo.t]* 


Here we recall d is the dimension of domain fh And we assume f 2* is a uniform 
rectangular grid in the above definition. In general the underlying discrete innerprod- 
uct must be modified, see Olsson [6]. In order to simplify the notation, a single index 
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j is used for the grid points in d dimensions. For the same reason, we use the same 
K and a in Definition 2 and 3. In fact, they are not necessarily the same constants. 
However, in order to approximate the original problem better, it is desirable to have 
the a in the discrete case to be close to the one in original problem. This is refered 
to as strictly stable in Olsson’s report. It can be shown that the two a’s are equal up 
to an error of order 0(h) in his estimates. 

Before deriving the error estimates for CAG methods, another assumption is nec- 
essary. We will assume K — 1, where K is the coefficient in equation (22). According 
to the stability results of Olsson, K > 1. Let’s see what happens if K > 1. Since 
adaptive grids are used, the stability bound in Definition 3 is used whenever a regrid- 
ding process is carried out. We consider a simple case with homogeneous boundary 
conditions and zero forcing term. Assume that we regrid at every time step, then we 
get 


IKWIIL < 

IK(2*)lln, < A-V“*IM0)||^ 


So at time 7 1 , the growth bound looks like 

IMr)||&. < tf TA e“ T IM0)||£„ 

This bound is useless unless K < 1 +0(k). However, we can always absorb K into the 
term e at with a new a if K < \ + 0(k). Therefore, we take K = 1. Fortunately, many 
problems we are interested in have this nice property since they describe physical 
phenomena which conserve energy. 


3.2 Convergence 

Now, we are ready to define in the case of CAG’s. We use the notation G to 
represent this discretization of the original domain. Given a composite component 
grid, which is defined in Section 2 and may contain several levels of subgrids, the 
domain where the computational solution vector Vh is defined can be constructed in 
the following way. Suppose there are L levels of subgrids Go, . . ., and Gl- i, and every 
component grid Gij is a regular grid. Then i ^ is defined spacially on the set 


G = Gi_i (J ( LJ (Gi-, - G,)) 


/=1 


(26) 
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where G/_i — Gi is the complement operation. This set is well defined in the 1-D 
case since we enforce Gi C G;_ i and there is no overlapping among the subgrids of 
the same level. But in higher dimensions there may be overlapping areas in Gi. The 
above definition is still valid if the solution vector is uniquely defined on each G;. 
All the components in G/ will have errors of the same magnitude. So the solution 
vector can be defined using any one of the grids in the overlapping area. Suppose 
G/jj, Gjj 2 , ..., Gij n have an overlapping area, where j x < j 2 < ... < j n - Then 
the solutions of G;,j a are used in the overlapping area. In other words, by using the 
smallest-indexed grid in the overlapping area, we can uniquely define the solution 
vector on every level. Then, equation (26) is well defined in higher dimensions and Vh 
is defined on a piecewise uniform mesh. We also notice that the set G is a function of 
time, i.e., G = G(t), because it may change whenever a regridding process is done at 
any level. By this observation, we can define Vh on the interval [0, T]. In fact, v k is a 
piecewise constant solution function on [0, T\. Its solution values change only when a 
integration or regridding process is carried out at any level. Because G is a piecewise 
uniform grid, we can define dG to be the union of all the boundaries of these uniform 
grids. The discrete norms defined in equations (23) to (25) can be easily generalized 
as follows: 

IM*)llo(«) = J2 ht u) \v h ( Xj ,t)\ 2 

x,eG(t) 

IIMIaG(T)x[0,T] = 

IIAIIg(T)x[0,T] = ^2 ^l(i)^Hj)\fh( x jiU)\ 

t,e [o,T] kXj eG(t) 

where h^j) is the mesh szie of the grid to which the grid point x } belongs. For 
simplicity, we assume the mesh sizes in all directions in a component grid are equal. 


Now we define the solution error e/* on the same domain of Vh as follow: 

efc(<) = v h (t) - u h {t) (30) 

where Uh is the projection of the exact solution onto G. 

On each section with uniform grid G um / orm and time interval [to,G], the error e* 
satisfies following error equations: 

(t T k ) = Lh. Gfc(^) T" k on G un if 0 rm ^ (31) 

Cfc(0) = Co/i G un ij OT mi (32) 

= ’b bdry O™ dG un if OTTn X [to,G]- (33) 


Since linear problems are considered here, the error equation is the same as the 
discrete equation for v h except the forcing and boundary terms are replaced by the 


(27) 

(28) 
(29) 
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local truncation errors. Thus, the growth bound for on the piecewise uniform grids 
can be obtained by applying the stability results in Section 3.1. 

Next we are going to prove the error estimation for CAG methods using induction 
on the number of grid levels. First we consider the case of only one level of grids. 
Since it is piecewise uniform, we can use the result in equation (22) on each piece 
of uniform grid. Also we bound the error on grid G by the summation of the errors 
on the pieces of uniform grid because of the way we define G and the norm || • ||g- 
Assume the local truncation errors of the numerical method used have the accuracy of 
order p at interior and boundary points. The initial values are also assumed accurate 
to order p. Then from equation (22), it is easy to see following error bound holds: 

IM^IIo < C„e» T C (34) 

where < means “asymptotically less than”, since only the leading terms of the trun- 
cation errors are estimated here. 

Now we assume that equation (34) holds for L levels of grids. Next we will show 
it for L + 1 levels. Because of the way we construct the subgrids which are nested 
in their parent grids, we only need to consider two levels. By induction, we know 
the error on the finer of the two levels is of order 0(h\ ), so we can always write it as 
O(fto). Assume the computation starts at t = we have the following error bound 
until we regrid at t = tj, 

| | e /i(^l)| |g ^ e0r ( tl ” t °)(l|e/i(^l)||G + I l^int | Igx [*i —t 0 ] + ll r Wry||aGx[ti-t 0 ]) (^5) 

where Ti nt is the local truncation error in the interior points, and r^dry is the truncation 
error introduced at the boundaries, which include the exterior boundaries dG exi and 
the interior boundaries dG xni . Suppose a pth order method is used in the interior 
points. On the exterior boundaries and the interior boundaries, suppose we use 9 th 
and rth order methods, respectively. Then the following estimates are obtained: 

llnn,|IL[.,- 1 „] £ C.*oV(G )((.-*<.), (36) 

lln*„III G)I | 11 - fcl < c,*W#<(0G“')(<i - <„) + c 3 kX r i‘(aa‘"‘)(t l - i„) (37) 

where /*(•) is a measure of the grid G and its boundaries which are defined as follows: 


p(G) 

- h t{j )’ 

x } eG 

(38) 

p(dG) 

= 

XjSdG 

(39) 


At t = ti, the regridding process is executed on this regular grid. Another type 
of error Tj ntt , which is due to the initialization of the subgrids using interpolation, is 
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introduced. Let’s assume this process has error of O(hl). Thus, after the regridding 
process at t = t\, we will have following error estimate: 

IM*i+C < IM*i)llo + lfc»X 

< \\e h (ti)\\ 2 G + C 4 h 2 0 s - (40) 

If equations (35), (36) and (37) are substituted into equation (40), and the fact that 
hi/ki = constant is used throughout computation, the total error after the regridding 
process at t = t\ is 

IMMIfe < + c,hly(G)(u - to) + c J A? ,+, ' f .(aG“ , )(t 1 - to) 

+C 3 A? r + ,, rfdGr‘)(t, - t„)f + C,hl‘. (41) 

where the subscript + on the left hand side of equations (40) and (41) represents the 
fact that the regridding process is done at time t = t\. All the terms in equation 
(41) are under control except fi(dG ,nt ) since we do not know how many subgrids are 
generated. The following lemma gives a bound on fi(dG >nt ). 

Lemma 1: Let Go be a regular grid in R d . Assume that there are L levels of grids 

G 0 , • • •, Gi-i with the refinement ratio m = hj/hj+\. Then the following asymptotic 
inequality holds for small h 0 . 

/*(UU dG u) £ h?C(d,m,L)r{G 0 ) (42) 

/=1 3 

where C(d,m,L) is a constant depending on d, m, and L. 

Proof: Since a regular grid is an union of rectangular grids, without loss of gen- 

erality, we can assume that Go is a rectangular grid. First we bound the boundaries 
of G\. According to the discussion in section 2, it can be shown that the case with 
largest interior boundaries as ho is very small for d = 2 is in Figure 4. So, it is easy 
to see that the following bound holds: 

MUS«u) S 2 dho' 1 ~ 

= v| 'e(Go). 

So, for each level /, l = 1,2, . . . , L — 1, we have 

m(U dGij) < hj\—^i(Gi-i) 

3 

, 2d . _ . 

^ h o oK^j) • 
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Figure 4: The case with largest interior boundaries in 2-D 


Therefore, we have 
L- 1 

U 

l=i j 




' 2 d ' 


• + (£) i_2 ) 


< h Q l C(d, m, L)h(Gq) 


where 


C(d, m, L) 


M(1 ~ 

2 J (1 - P) ' 


(43) 


Remark: The above lemma tells us is that, under the worst senerio, the error term 

involving fi(dG ,nt ) may lose one order of accuracy. However, the bound in the lemma 
is very pessimistic. In most situations, the number of subgrids is usually quite small. 
So it is very unlikely the term n(dG xnt ) will actually change the order of accuracy at 
the interior boundaries. Also we see that it is not a good idea to use too many levels 
with large mesh ratio, since the constant C in Lemma 1 will be very large when L is 
large and p > 1. 

With this lemma, we can finally bound the error. We notice that for our current 
proof, we only need the bound for two level case, i.e., 

MU aG lj) s (44) 

3 

Suppose the worst case is considered here, i.e., the refinement and reinitialization 
processes are assumed as often as possible. Assuming the time interval is [0, T], then 
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we have following error bound using equation (41) and (44): 

IKCnilc S e" T (IWo)||J + c I CrMG) + c 2 C ,+,) T#.(aG'«) 

+C 3 k?’ +,| r^V(Go) + (45) 

\{ q = p — 1 ? r = p, s — p and the initial error is order 0(h q), then equation (45) can 
be written as 

IMr)||c < Ce^K (46) 

where the constant C may depand on T, the number of levels. Therefore, by the 
induction argument, we have proved the above error bound for CAD’s. Convergence 
follows as h 0 0 with a fixed maximum number of levels. We have just proved the 
following theorem. 

Theorem 1: Suppose a well-posed initial-boundary value problem (equations (1) - 

(3)) is solved by a numerical method (equations (4) - (6)) using component adaptive 
grids with a fixed maximum number of levels of grids. Also assume that the stability 
condition (equation (22)) for the numerical method holds on uniform grids, and the 
accuracies for the numerical approximations are order p, p— 1, and p on interior points, 
exterior boundaries and interior boundaries, respectively. Also that the interpolation 
process used in the CAG method has accuracy of at least order p. Then the numerical 
solution converges to the exact solution in the norm || * \ \g as h 0 — j ► 0, and the order 
of convergence is p. 


3*3 Error Estimation 

From the discussion in Section 3.2, the sources and the magnitudes of various com- 
putational errors are well understood. We can implement the numerical method so 
that the local truncation error r tn * in the interior points is dominant. Then 

|Mr)||c < \A y, 'V(G) T max(r jnl (j, i)) 

< £. (47) 

Therefore, if the bound 8 for the local truncational error satisfies following relation 

6 = . e , (48) 

y/^MG) T 

then the final error is guaranteed to be bounded by e. 


max(r; nt (x,<)) < 
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Not only can we bound the error, we can also estimate the final error. Since local 
truncation errors in all grid points are estimated once every several time steps, simple 
quadrature can be used to estimate the final error. In many applications the growth 
factor a is not known in advance but it is very easy to approximate a using the 
computed solution as the growth factor is asymptotically the same for the solution 
and the error. 


4 Computational Examples 

In this section some numerical experiments in one and two space dimensions are pre- 
sented to illustrate component adaptive grid methods and the error bounds derived 
in section 3. The programs are written in C because of its capability for dynamic 
memory allocation and flexible data structures, which are crucial for our CAG meth- 
ods. A SUN SparclO workstation with 96 Mbytes memory is used for both the 1-D 
and 2-D cases. 

Example 1 (1-D wave equation). In this example, we compute the solution to 
the following wave equation in one space dimension. 


u t 

u(x, 0) 


— u x + cm x € [0,1], t e [0,0.4] 

I 0 if X€ [0,0.2] U [0.4,1] 

( f(x) if a: € (0.2, 0.4) 


where 

/(*) = 50 (0.01 - (x - 0 .3) 2 exp(^ 7 -^ML-_ T ). 

The exact solution is a wave front traveling from left to right with speed 1 and growth 
rate e at , where a is a constant. Two methods are used to solve this equation: the first 
order up-wind method and second-order Lax-Wendroff method. Second order Hermite 
interpolation is used with both methods. In Figure 5, we plot the exact solution and 
the computational solutions at t = 0.4 with a = 0 using the Lax-Wendroff method 
on the coarse grid, 3-level adaptive grids with mesh ratio m = 4 and uniform grid 
with mesh size equal to the smallest of the adaptive grids. Both the solutions on the 
fine and 3-level adaptive grid are much better than the one on the coarse grid, which 
has wiggles at the left corner where the local truncation errors are large. However, 
the adaptive grid uses much less time than the uniform fine grid. 

Next we use the results in section 3 to control the local truncation error tolerance 
8 according to the final error bound e. All the computations here are done with 
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0.2 


exact solution 


coarse solution (h=0.01) 



3 level adaptive solution 




fine solution (h=0.01/16) 



Figure 5: various computional results to 1-D wave equation at t=0.4 


following parameters. 


mesh ratio 

4 

buffer zone width 

4 points 

regrid every 

16 steps 

coarse mesh 

0.01 

CFL No. A 

0.9 

growth factor a 

0 

final time T 

0.36 


We collected the data in Tables 1 and 2. In the first column, the final error bounds 
are given. The number of levels of adaptive grids used during computation is listed in 
the second column. The exact error and estimated error using simple quadrature are 
listed in columns 3 and 4, respectively. In the last two columns, we put the running 
times for adaptive grids and uniform grids with mesh size equal to the finest mesh 
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size in the corresponding adaptive grids. 


Table 1 


Results using the up-wind method for the 1-D wave equation 


6 

levels 

exact ||eft||c 

est. ||e fc || G 

time 

(sec) 

time using 
fine grid (sec) 

5 x 10~ 3 

i 

4.57 x 10" 3 

4.54 x HT 3 

0.0 

0.0 

1 x HT 3 

3 

9.37 x HT 4 

9.26 x 10- 4 

0.4 

4.7 

5 x 10~ 4 

3 

3.88 x HT 4 

3.82 x 1(T 4 

1.0 

4.7 

1 x 10" 4 

4 

9.98 x 1(T 5 

9.90 x 10~ 5 

12.9 

77.0 


Table 2 


Results using the Lax-WendrofF method for the 1-D wave equation 


e 

levels 

exact He/Jlc 

est. He/Jlc 

time 

(sec) 

time using 
fine grid (sec) 

1 x 10- 3 

2 

2.08 x 10" 4 

1.88 x 10~ 4 

0.1 

0.4 

5 x 10~ 4 

2 

1.86 x 10“ 4 

1.83 x 10" 4 

0.1 

0.4 

1 x 10" 4 

3 

3.68 x 10“ 5 

3.81 x 10" 5 

0.8 

7.0 

5 x 10" 5 

3 

2.24 x 10" 5 

2.15 x 10" 5 

1.0 

7.0 

1 x 10~ 5 

4 

7.75 x 10- 6 

7.70 x 10“ 6 

4.6 

115.0 

5 x 10- 6 

4 

2.68 x 10" 6 

2.60 x 10" 6 

10.0 

115.0 


Several interesting facts are illustrated in Tables 1 and 2. First of all, we see that 
our adaptive strategy is very efficient for solving PDE’s. It does efficiently generate 
different subgrids in response to the final error tolerances. For example, when we use 
Lax-WendrofF with tolerance e = 1 x 10 -3 and e = 5 x 10 -4 , two levels of grids are 
used in both cases. However, in order to satisfy the final error tolerance, the two 
Gi’s are constructed differently. This is shown in Figures 6 and 7. In the case of 
c = 1 x 10~ 3 , two small subgrids are generated around the two corners where large 
local truncation errors appear. When e is reduced to 5 x 10 -4 , a large subgrid is 
created in G\. The running times in Tables 1 and 2 show that the speedup, i.e., the 
ratio between the time using a uniform fine mesh and the time using an adaptive 
grid, increases as the final error tolerance is decreased. In other words, our adaptive 
strategy is more attractive when high accuracy is needed. This is because only very 
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small regions (the two corners in this example) need to be refined. The data also 
illustrate that our simple quadrature error estimation formula gives very satisfactory 
results. Of cause, one reason we have such accurate estimates is that the equation has 
constant coefficients. A more realistic variable coefficient problem is considered in the 
next example. As mentioned in the Section 3, the errors of interpolation and the ones 
at boundaries are assumed to be relatively small compared to the local truncation 
errors. Now we compute these two types of error explicitly and list them in Tables 
3 and 4. To be consistent with the notation used in Section 3, we use | |e tnt *t | |g? to 
represent the subgrids’ initialization errors caused by interpolation of the coarse grids’ 
data. | l^int | |g and || e 6dry||aG are used to represent the errors caused by local trucation 
errors on the interior points, and the errors on both exterior boundaries dG ext and 
interior boundaries dG tn \ respectively. Indeed, it is shown that | (e^nttl |c? and ||e^ ry ||aG 
are negligible compared to ||e m t||G- Finally we look at error estimation for different 
o’ s, since for some nonlinear problems, the error equations contains non-differentiable 
terms. The results for different a’s using Lax-Wendroff with tolerance e = 0.0005 are 
shown in Table 5. The estimated growth factors are listed in the table. Our error 
control and estimation works very well for all three test cases, a = 0, 3 and 6. 


G u Gi i2 

i i i 

Gcu 


Figure 6: Adaptive grid structure for e = 0.001 


j 


Gi,x 


G o,i 

i 

Figure 7: Adaptive grid structure for c = 0.0005 
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Table 3 


Various errors using the up-wind method for the 1-D wave equation 


e 

levels 

exact | e; nt ||G 

est. ||6|n*||c? 

€St. 1 1 Zbdry \ \dG 

6St. ||C|' nt 't||G 

5 x lO’ 3 

1 

4.57 x 10" 3 

4.54 x 10" 3 

5.25 x 10" 6 

2.13 x lO" 11 

1 x 10" 3 

3 

9.37 x 10" 4 

9.26 x 10~ 4 

5.55 x 10"’ 

3.27 x 10" 6 

5 x 10" 4 

3 

3.88 x lO" 4 

3.82 x 10~ 4 

7.96 x lO" 7 

9.16 x lO" 7 

'rr 

1 

o 

1 — 1 

X 

4 

9.98 x 10" 5 

9.90 x 10“ 5 

6.89 x 10~ 8 

2.41 x lO" 7 


Table 4 


Various errors using the Lax-Wendroff method to the 1-D wave equation 


e 

levels 

exact 

M 

G 

est. | |e tn t | |c; 

est. ||e(,drj(||aG 

est. | 

Cinit G 

1 x 10“ 3 

2 

2.08 x 10" 

4 

1.88 x 10" 4 

1.01 x 10" 6 

1.92 x 10" 7 

5 x 10~ 4 

2 

1.86 x 10" 

4 

1.83 x 10" 4 

1.04 x 10" 9 

7.80 x lO” 12 

1 x 10" 4 

3 

3.68 x 10“ 

-5 

3.81 x 10" 5 

1.43 x 10- 7 

3.17 x 10- r 

5 x 10" 5 

3 

2.24 x 

10" 


2.15 x 10 -5 

7.08 x 10“ 8 

8.03 

X 10" 8 

1 x 10“ 5 

4 

7.75 x 10" 

IT - 

7.70 x 10" 6 

1.31 x 10“ 8 

1.01 

x 10~ 7 

5 x 10" 6 

4 

2.68 x 10" 


2.60 x 10" 6 

4.65 x 10" d 

4.06 x 10" 8 


Table 5 

Estimates using the Lax-Wendroff method 
with various o’s and t — 5 x 10 -4 


a 

estimated a 

levels 

exact lle^llG 

est. He^llG 

0 

-1.43 x 10~ 4 

2 

1.86 x 10" 4 

1.83 x 10“ 4 

3 

2 + 2.56 x 10" 4 

3 

2.25 x 10" 4 

2.35 x 10" 4 

6 

6 + 4.39 x 10" 4 

4 

1.05 x 10" 4 

1.23 x 10" 4 


Example 2 (2-D rotating cone). The rotating cone problem has been used by 
Berger and Oliger [2] to illustrate the adaptive grid method using rotated rectangular 
subgrids. The problem is 
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where 


u t = yu x — xu y x € [— 1, 1], y € [—1, 1] and t € [0, 3.125] 

JO if (x-0.4) 2 + 1.5y 2 >0.04 

°) “ { f(x) if (x -0.4) 2 + 1.5y 2 < 0.04 

f{x) = 0.25 (1 - 25((x - 0.4) 2 + y 2 )) 2 . 

The solution is a cone with elliptical base which rotates counterclockwise about the 
origin. The solution is integrated until the cone is approximately halfway through the 
first revolution. The Lax-Wendroff method and second order interpolation are used in 
this example. The boundary conditions are zero inflow and first order extrapolation 
at outflow. All computations are done with following parameters: 


mesh ratio 

4 

buffer zone width 

2 points 

regrid every 

8 steps 

coarse dx 

0.05 

coarse dy 

0.05 

CFL No. A 

0.25 

final time T 

3.125 


Snapshot views of the one subgrid at three intermediate time steps are shown in 
Figure 8. We should mention that this example does not fully explore the power of 
our stair step subgrids, since the refined regions here are elliptical, which are easily 
covered by rectangular grids. We expect our stair step subgrids to be more efficient 
when the geometries of subgrids are more complicated. Various computational results 
are plotted in Figure 9. It shows that the solution using a coarse grid has lots of 
oscillations. They are reduced dramatically if 2 levels of adaptive grids are used. If 
we use 3 levels, the oscillations are invisible. Also, we see that the solutions using 
adaptive grids are as good as those using uniform grids with the smallest mesh size of 
the corresponding adaptive grid, which is shown by their errors (see Table 6). Next, 
our error estimation developed in Section 3 is used to approximate various types of 
computational errors. The results are listed in Tables 6 and 7. 

Since we are solving a variable coefficient problem, the estimated ||cfc| |c? is slightly 
larger than the exact one. Also, like the 1-D wave equation example, the local trun- 
cation error 1 1 e^ nt 1 1 g dominates the other two types of error, \\e b dry\\dG and \\e in i t \\G 
(see Table 7). Finally, we estimate the efficiency of our stair step subgrids. We use 
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the two level case as an example. The area of the subgrid is about 16% of the original 
domain. So, the running time on this subgrid is about 167 x 16% = 26.7 (sec). Also, 
we know the running time on the coarse grid is 2.7 (sec). Thus, the total time spent 
on other things besides integrating the grids is about 31.7 — 26.7 — 2.7 = 2.3 (sec), 
which is approximately 8% of the total time for this two level adaptive computation. 
This percentage is less than that reported in [2], which was about 12% when using 
rotated rectangular subgrids. Like our 1-D example, the 2-D example also shows that 
the speedup increases as we decrease the final error tolerance. As mentioned earlier, 
we expect our stair step subgrids to be more efficient when refined regions with more 
complicated geometries are encountered. 



x 


Figure 8: Stair step subgrids for the rotating cone problem 
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exact 



i T 0 _1 


X 


2 level adaptive 



3 level adaptive 



Figure 9: 


Solutions for 


coarse (dx=dy=0.05) 



uniform (dx=dy=0.01 25) 



uniform (dx=dy=0.003125) 



cone problem 



Table 6 


Results using the Lax-Wendroff method for the 2-D rotating cone problem 


e levels 

exact ||efc||G 

est. ||e fc || G 

time 

(sec) 

time using 
fine grid (sec) 

exact ||e /l || G 
for fine grid 

■DEBI 

2 

5.11 x 10~ 3 

6.70 x HT 3 

31.7 

167.0 

5.11 x 10“ 3 

i x nr 3 

3 

5.32 x 10" 4 

7.40 x 10" 4 

1024.0 

10879.0 

5.27 x IQ’ 4 


Table 7 

Various errors using the Lax-Wendroff method for the 2-D rotating cone problem 


e 

levels 

exact | eh | G 

est. ||e int || G 

est. ||efcd ry ||a G 

est. ||e,' n ,‘t||G 

m 

6 

2 

5.11 x 10~ 3 

6.70 x 10" 3 

4.77 x 10~ 6 

1.94 x 10" 4 

1 x 10- 3 

3 

5.32 x lO" 4 

7.40 x 10" 4 

2.96 x 10-* 

2.06 x lO" 5 


5 Conclusions 

We have presented an algorithm for solving PDE’s and estimating the final error 
with component adaptive grid methods. Good efficiency of such grids has been illus- 
trated in our 2-D example. Applying recent stability results [6] for a class of numerical 
methods on uniform grids, we have proven the convergence of these methods for linear 
problems on CAG’s. We obtain satisfactory final error estimates. Since local trunca- 
tion error estimates are obtained during the process of grid refinement, the amount 
of extra work to obtain a final error estimate is very small. For linear problems, not 
only can the algorithm estimate the final error, it can also be used to control the tol- 
erance for the local truncation error in terms of the given final error bound. The next 
challenge is to estimate errors for nonlinear problems. Although it is very hard to 
establish such results for general nonlinear PDE’s, Olsson and Oliger [7] have devised 
a technique that makes it possible to obtain energy estimates for initial-boundary 
value problems for a class of nonlinear conservation laws. The estimates for finite 
difference methods on such nonlinear problems is currently under development. We 
think such results will help us estimate and control errors on our CAG’s for a large 
class of nonlinear problems. 
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